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PREFACE 

This  Memorandum  is  part  of  RAND's  continuing  search  for  new 
ways  of  utilizing  the  modern  digital  computer.  The  authors  present 
a  method  for  numerically  integrating  nonlinear  partial  differential 
Integral  equations,  which  occur  in  such  fields  as  radiative  transfer 
and  mathematical  biology.  The  method  is  then  specifically  applied 
to  solving  a  basic  equation  of  transport  in  a  spherical  shell. 
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SUMMARY 

In  this  Memorandum,  the  authors  show  how  to  approximate  a  non¬ 
linear  partial  differential  integral  equation  by  a  system  of  ordinary 
differential  equations.  A  table  of  necessary  constants  is  provided, 
and  the  results  of  a  test  calculation  on  an  equation  of  radiative 
transfer  in  a  spherical  shell  are  described. 
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I.  INTRODUCTION 


In  applying  invariant  imbedding  to  the  radiative  transfer  proc 
esses  associated  with  plane  parallel  regions,  we  encounter  a  func¬ 
tional  equation  of  the  form 


||(t,v,u)  +  (^  +  £)s  -  g(u,v,*)  (1) 


1 

1 

+  X 

l+-ij 

.  ,  1  r  a/  •  \ 

1  +  2  J  S(*,v  *u)  “T7 

o 

o 

- 

- 

where  S(o,v,u)  ■  0  (see  Refs.  1,  2,  and  3).  Here  z  2  0,  0  i  u,  v  J  1. 
This  can  be  approximated  by  means  of  a  finite  dimensional  set  of 
ordinary  differential  equations  by  introducing  quadrature  techniques. 
Write 


1  N 

J  S(s  ,v,u ')  —t  ^  £  w1S(r,v,xi)/xi  ,  (2) 

o  i-1 

1  N 

J  S(*,v',u)  “"“7  a  £  w1S(s,x1,u)/x1  , 

o  i-1 

where  XpX^ . x^  are  ^  roots  of  the  shifted  Legendre  polynominal, 

P„(x)  -  P„(l-2x).  Using  these  approximate  relations  and  setting 
S(z,x^,x  )  ■  S^(x),  Eq .  (1)  reduces  to  a  finite  dimensional  system 
subject  to  initial  conditions.  This  technique  has  be  quite  success¬ 
ful  in  practice,  as  evidenced  by  the  results  in  the  cited  references. 
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If  we  turn  to  the  study  of  corresponding  trsnsfer  processes  for 
spherical  and  cylindrical  regions,  we  meet  a  much  more  formidable 
equation 


as,  N 

5l(z’v-u> 


1— u 
zu 


2  2 


as  .  n  .  i  »  „  +u  \  s 

st +  Cv +  ~J s  -  K~n)  ~z 


v  u 


(3) 


1 

1 

■  g(u,v,z)  +  \ 

1  +  |  J  S(z.v.u')  ij. 

1  +  r  S(Z,V  ',u) 

^  si  V 

o 

0 

In  the  following  we  shall  briefly  sketch  an  approximation  technique 
which  enables  us  to  reduce  the  numerical  solution  of  Cq .  (3)  to  that 
of  a  finite  system  of  ordinary  differential  equations,  and  shall  also 
describe  a  sample  calculation.  Mora  detailed  results  will  be  presented 
subsequently.  The  method  has  been  applied  successfully  to  a  number 
of  other  classes  of  functional  equations  involving  partial  derivatives. 
Finally,  we  note  that  equations  involving  partial  derivatives  ar d 
integrals  occur  with  great  frequency  in  mathematical  biology.^ 

An  approach  of  Chandrasekhar’s  is  described  in  Ref.  1. 
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II.  DERIVATIVES  AS  LINEAR  COMBINATIONS  OF  FUNCTIONAL  VALUES 

In  order  to  generalize  the  approach  to  Eq.  (1),  we  replace  the 

partial  derivatives  S  and  S  by  linear  combinations  of  the  values 

u  v 

of  S  at  the  points  u,v  ■  x^,  x^,  i,j  ■  1,2,...,N.  Given  a  function 
f(x),  we  want  an  approximate  relation  of  the  form 

N 

f/(xi>  -  £  atJ  f(x  )  ,  i  -  1,2 . N  (4) 

j-1 


To  determine  the  coefficients  a^ ,  we  ask,  by  analogy  with  the  Gaussian 

quadrature  formula,  that  Eq.  (4)  be  exact  if  f (x)  is  a  polynominal 

of  degree  N  1  or  less.  To  obtain  a^j,  we  use  the  test  functions 

f  (x)  -  PM(x)/[(x-x  )  P  '(x  )1.  A  simple  calculation  then  yields 
m  N  L  m  N  m  J 


pn'<V 


ira 


<VXm>  PN  (xm> 


i  f  in 


(5) 


PN  <XB> 


(l-2x  ) 
m 


nxn 


2P  '(x  )  2(x2-x  ) 

N  m  'mm' 


for  m  ■  1,2,...,N.  In  view  of  the  symmetry  of  P  (x)  about  x  ■  .5,  it 

N 


is  clear  that  -  -a 


N+l-J’  a  re,u^t  yields  both  a 

useful  check  on  the  calculation  of  these  parameters  and  a  reduction 


in  the  size  of  the  tables.  A  table  of  values  of  a^  for  N  ■  3,7,9 
follows.  These  values  were  calculated  by  H.  Kagiwada  and  verified 


by  J.  Jolissant. 
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Table  1 

THE  COEFFICIENTS  FOR  N  *  5 


i  -  1 

-0. 10134081E  02 
-0. 11035337E  01 
1-2 

-0. 19205120E  01 
0. 48883323E-00 
i  -  3 

0.60233632E  00 


0. 15403904E  02 

-0. 15167064E  01 

-0. 2870 7765E  01 


-0. 80870874E  01 

0.48055013E  01 

-0. 35527137E- 14 


0. 39207982E  01 

-0. 18571160E  01 

0. 28707765E  01 


Table  2 

THE  COEFFICIENTS  a  FOR  N  -  7 

i  -  1 


-0. 19 136364E 

02 

0. 30166068E 

02  -0. 18343136E 

02 

0. 12020668E 

02 

-0. 735540 34E 

01 

0. 37037V09E 

01  -0. 10536210E 

01 

1  *  L 

-0.30  7  7400  IE 

01 

-0. 32947313E 

01  0. 94826608E 

01 

-0.49 141 384E 

01 

0. 27743267E 

i  _  *1 

01 

-0. 1348560VE 

01  0. 37784329E- 

GO 

1  "  J 

0.  738786V IE 

00 

-0. 37433740E 

01  -0.971  74703E 

00 

0. 56413488E 

01 

-0. 24639939E 

01 

0. 10V51929E 

Cl  -0. 29621352E- 

00 

1  ■  4 

-0. 36940283E- 

■00 

0.  14803137E 

01  -0. 43048331E 

01 

-0.99475983E- 

13 

0. 4304  8  331E 

01 

-0. 14803137E 

01  0. 36940283E- 

00 
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Table  3 

THE  COEFFICIENTS  a^  FOR  N  -  9 


i  -  1 

-0. 30899183E  02 
-0.1 6634 32 5E  02 
-0. 10328869E  01 
1  -  2 

-0. 4632184 7E  01 
0. 589 5008 7E  01 
0  .  33923594E-00 
1  ■  3 

0. 9977960 8E  00 
-0.46474057E  01 
-0. 22383800E-00 
i  -  4 

-0. 41927865E-00 
0. 67044574E  01 
0. 20889 316E-00 
i  -  5 

0. 25654 308E-00 
0. 56843419E- 11 
-0. 25654308E-00 


0. 49462 602E  02 
0. 11463908E  02 


-0.  5554064 7E  01 
-0.  39077266E  01 


-0. 51953604E  01 
0. 27969636E  01 


0. 172 38 12  3E  01 
-0.  308400  7 5E  01 


-0. 9  70  60200E  00 
0.  56744949E  01 


-0. 31847722E  02 
-0. 71444762E  01 


0. 15529632E  02 
0. 23856884E  01 


-0. 1966641 7E  01 
-0.  1630  3  3  3  5E  01 


0.  23009713E  02 
0.  362 2 3 7 1  IE  01 


0.88594615E  01 
0.  1196127 7E  01 


0. 90706996E  01 
0. 7981200 6E  00 


■0. 5275564 3E  01  -0.  72470224E  00 
0. 16267280E  01  -0. 76033?20E  00 


0. 22877170E  01  -0.56744949E  01 
■0.2287 7 l 70E  01  0.97080200E  00 
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III.  SAMPLE  CALCULATION 

Presented  below  are  the  results  of  a  particular  calculation  of 
interest  in  determining  the  flux  reflected  by  a  spherical  shell 
atmosphere.  We  approximated  Eq.  (3),  with  g(u,v,z)  =  0,  by  the 
system  of  ordinary  differential  equations 


with  initial  conditions  S^j(a)  *  0. 

We  set  \  -  1,  usually  the  severest  test,  and  integrated  over  z 
from  a  to  a  +  3  with  the  values  a  -  100,  500  and  1,000.  The  constant 
a  is  the  inner  radius  of  the  spherical  shell.  For  comparison  with 
the  plane  parallel  case  in  Ref.  1,  we  also  printed  out  values  of 
rij(z)  -  SiJ(z)/Ax1. 

Some  results  are  shown  graphically  in  Fig.  1.  Note  that  as 
the  radius  of  the  inner  surface  of  the  shell  increases,  the  reflection 
function  of  the  shell  approaches  that  of  the  slab.  This  is  one  test 


Angle  of  reflection 


Fig.l  —  Some  reflected  intensity  patterns  for  shells  with  albedo 
X  =  1  and  thickness  x  -  3,  for  various  angles  of  incidence 
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of  the  validity  of  the  approximation  technique  used.  In  addition, 
comparison  of  the  N«7  and  N*9  cases  indicates  excellent  agreement. 
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